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ABSTRACT 

Recent multi-scale simulations have made it possible to follow gas inflows responsible 
for high-Eddington ratio accretion onto massive black holes (BHs) from galactic scales to 
the BH accretion disk. When sufficient gas is driven towards a BH, gravitational instabilities 
generically form lopsided, eccentric disks that propagate inwards from larger radii. The lop- 
sided stellar disk exerts a strong torque on the gas, driving inflows that fuel the growth of the 
BH. Here, we investigate the possibility that the same disk, in its gas-rich phase, is the putative 
"torus" invoked to explain obscured active galactic nuclei and the cosmic X-ray background. 
The disk is generically thick and has characteristic ~ 1 — 10 pc sizes and masses resembling 
those required of the torus. Interestingly, the scale heights and obscured fractions of the pre- 
dicted torii are substantial even in the absence of strong stellar feedback providing the vertical 
support. Rather, they can be maintained by strong bending modes and warps/twists excited 
by the inflow-generating instabilities. A number of other observed properties commonly at- 
tributed to "feedback" processes may in fact be explained entirely by dynamical, gravitational 
effects: the lack of alignment between torus and host galaxy, correlations between local SFR 
and turbulent gas velocities, and the dependence of obscured fractions on AGN luminosity 
or SFR. We compare the predicted torus properties with observations of gas surface density 
profiles, kinematics, scale heights, and SFR densities in AGN nuclei, and find that they are 
consistent in all cases. We argue that it is not possible to reproduce these observations and the 
observed column density distribution without a clumpy gas distribution, but allowing for sim- 
ple clumping on small scales the predicted column density distribution is in good agreement 
with observations from Nu ^ 10 20 — 10 27 cm -2 . We examine how the Nu distribution scales 
with galaxy and AGN properties. The dependence is generally simple, but AGN feedback 
may be necessary to explain certain trends in obscured fraction with luminosity and/or red- 
shift. In our paradigm, the torus is not merely a bystander or passive fuel source for accretion, 
but is itself the mechanism driving accretion. Its generic properties are not coincidence, but 
requirements for efficient accretion. 
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1 INTRODUCTION 

It has long been realized that bright, high-Eddington ratio accretion 
(i.e., a quasar) dominates the accumulation of mass in the super- 
massive BH population (Soltan 1982, Salucci et al. 1999, Shankar 
|et al.| 2004 , Hopkins et al. 2006d ). The discovery, in the past decade, 
of tight correlations between black hole mass and host spheroid 
properties including mass ( Kormendy & Richstone 1995 , Magor- 
|rian et al.]|1998| ), velocity dispersion (Ferrarese & Merritt 2000, 
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Gebhardt et al. 2000 ), and binding energy or potential well depth 
([Hopkins et al.|2007b|a| |Aller & Richstone|2007||Feoli & Mancini| 
2009 ) implies that black hole (BH) growth is tightly coupled to the 
process of galaxy and bulge formation. Increasingly, models invoke 
feedback processes from active galactic nuclei (AGN) to explain a 
host of phenomena, from the origin of the M B h — & relation, to rapid 
quenching of star formation in bulges, to the buildup of the color- 
magnitude relation and resolution of the cooling flow problem (see 



e.g.|Silk&Ree s 1998 , P i Matteo etaT]2 005 , Springel et al. 2005a 



Hop kins et al.|2008al [Hopkins & Elvis |2010HCroton et al.|2 006 
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Cattaneo et al. 2009 and references therein). 
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Observations have demonstrated that most of the accretion lu- 
minosity in the Universe is obscured by large columns of gas and 
dust (e.g.[Lawrence||199l||Risaliti et al.|[T999] [ffill et al.|[1996 l 
Sim pson et al.|1999[|Willott et al.|20001|Simpson & Rawlings 2000 
Hao et al. 2005 ; Ueda et al.|2003[ and references therein). This ob 
scured AGN population dominates the population of X-ray sources 
(|Miyaji et al.|2001[|Ueda et al. 12003] |Nandra et al. 12005 [|Hasingi 
et al.|2005[|Steffen et d.|2003||Grimes eTa?T2004, Hasinger 2004 



Sazonov & Revnivtsev||2004||Barger & Cowie||2005[ |Gilli et at] 



2007 ; Hasinger 2008), and accounts for most of the observed X-ray 
background luminos ity ([Setti & Wo ltjer 1989 ||Madau et al.|1994l 
[Comastri et d.|1995||Treister & Urry|20061|Gilli et al.|2007| ). It may 
dominate the bright end of the infrared luminosity function as well 
([Sanders & Mirabel|1996[|Komossa et al.|2003[|Ptak et al.|2003 



Hicko x et al.|2007||Daddi et al.|2007 [[Alexander et aL|2 008 ; Hop- 
kins & Hernquist 20101). The abundance of obscured quasars re- 
mains a major uncertainty in reconciling synthesis models of AGN 
populations with the BH mass function today, and (by implication) 
understanding the radiative efficiencies of quasars ([Salucci et al. 



1999 



2009 



|Yu & Tremaine|2002[ [Hopkins et al.|2007c[[Shankar et al. 



Various specific galaxy populations (for example EXOs, 



XBONGs, ULIRGs, SMGs) commonly host obscured AGN {[Yuan 



& Narayan|2004[[Geo rgantopo ulos & Geo rgakakis 2005 , |Max et al. 



2005 



2007 



2009 



Maini eri et al.||2005] [Alexander et al 
Alexander et al.||2008 



,2005, Daddi et al. 
Riechers et al.J[2008[ |Trump et al. 



Georgakaki s et al.||20Q9| |Nardini et aL] [2009 ). And theo 



retical models have long predicted that in violent events such as 
galaxy mergers, there should be a transition from an early, "buried" 
accretion stage corresponding to e.g. "warm" ULIRGs and sim- 
ilar galaxies, to an at least partially un-obscured phase in which 
the AGN removes dust and gas and is visible as a bright quasar 
([Sanders et al.|1988a|b[[King|2003[|Di Matteo et al.|2005l|Hopkins| 



etal.|2006a|c|b||200 8b 2006e| |2010[|Hopkins & Hernquist| 2010 
Younger et al.|2009| >. 

Yet AGN obscuration remains poorly understood. The most 
popular models invoke a torus-shaped "donut" of obscuring mate- 
rial, on scales anywhere from ^ 0.1 — lOOpc, to explain most of 
the heavily obscured AGN population ( Antonucci 1993 ; Lawrence 
|1991| >. If one empirically assumes unification of obscured and un- 
obscured AGN, then a number of the properties of the torus can 
be inferred: scale radii somewhere in the range above, and scale 
heights h/R of order ~ 1/3 (Risaliti et al. 1999). The observed dis- 
tributions of quasar and AGN column densities, and their detailed 
SED properties, place strong constraints on the densities, structure, 
and column densities within the obscuring material, with typical 
column densities as high as ~ 10 26 cm -2 through the edge-on plane 
of the material. And direct observations are beginning to probe 
these scales, through combinations of diverse techniques such as 



adaptive optics and maser observations ( Greenhill et al. 2003 , Jaffe 



|et~aT 2004, M ason et al. 2006, Sanchez^^j2006, Davies et al. 
[20061 IKrips ct al.||2007l |Davies et al.|[2007l |Hicks et al.|[2009| 
Ramos Almeida et al. 2009), giving constraints on the kinematics, 
gas and dynamical masses, and star formation rates at these scales. 
Indeed, this simple model of obscuration has proven successful at 
explaining a large number of AGN observables, and the torus forms 
the basis of most models uniting Type 1 and Type 2 AGN. 

These successes should not mask the fact that the torus re- 
mains a phenomenological model. The simple "donut" picture is 
just a toy model - there are a large and growing number of un- 
ambiguous cases where it fails, whether in predicting detailed 
radiative transfer properties coming from the microphysical gas 
structure (Mason et al. 2006, Elitzur & Shlosman 2006, Mor et al. 



)6l J20( 
)0| |20( 



;er 
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2009), or where the implied torus properties would involve bizarre 
radii and/or dust temperatures ( Kuraszkiewicz et aT]|2000[ |Tran 



2003l|Page et al.|2004[ [Stevens et al.|2005[|Ramos Almeida et~ 
2009), or where it is simply clear that the dominant obscuration 
is isotropic, or time dependent, or comes from much larger scales 
(e.g. those associated with circumnuclear starbursts and/or the host 
galaxy; see|Soifer et al.|1984{ |Scoville et al.|1986[ |Sanders et aL 



1988a b , Zakamska et al. 2006 ;|Liu et al.|2009[|Donley et al.|2 005 , 
Rigb y et al.|2006[|Sc hartman n et al.|2005[|Hatziminaoglou et al.| 



2009, Rowan-Robin sonet al.|2009[|Martinez-Sansigre et al.|2009[ 
Lagos etal.|2011| ). 

Without a physical model for the origin and evolution of nu- 
clear gas inflows, a large number of questions remain unanswered. 
Where do toroidal-like obscuring gas structures come from, in the 
first place? What determines their characteristic gas masses, radii, 
and structure? Why are such structures ubiquitous around AGN? 
Are they, in fact? It is also usually assumed that the torus is simply 
an obscuring "bystander" to the accretion event, or at most a passive 
fuel reservoir. But could the torus play some more critical role in 
the accretion process itself? A major long-standing puzzle is what 
drives and maintains the scale height of the torus - simple thermal 
pressure would be lost to cooling in a time much shorter than the 
local dynamical time. A large number of torus properties have been 
attributed to feedback from either young stars or the BH accretion 
itself - including the typical scale heights, clumping/phase struc- 
ture, gaseous velocity dispersions, possible correlations between 
these quantities and star formation, and even the fueling rates onto 
the BH (e.g. |Wada & Norman|2002l |Schartmann et al.|2009[ and 
references therein). But it is important to recall that we do not yet 
understand the basic dynamics of gas and stars entirely in the ab- 
sence of feedback! 

There have been some attempts to address these from a phys- 
ically motivated perspective, both in analytic and numerical work 
([Kawakatu & Wada|2008[|Cattaneo et al.|20051|Elvis|2000[|Hop^] 
|kins & Elvis|2010[|Elitzur & Shlosman|2006[ |Wada et al.|2009] L 
However, analytic models are severely limited by the fact that the 
systems at these radii are highly non-linear, often chaotic, and not 
necessarily in steady state (with inflow, mass buildup, star forma- 
tion, and feedback processes all competing). If one wishes to simul- 
taneously follow the torus itself and the chaotic, non- symmetric gas 
inflows that form it in the first place, simulations are necessary. But 
simulations of galactic scales used to follow inflows and AGN ob- 
scuration have resolution of ~ lOOpc, much larger than the relevant 
scales here (Cattaneo et al. 2005 , Hopkins et al. 2005a b). 

Although progress has been made with zoom-in refinement 
techniques (see e.g. Escala 2007 , Levine et al. 2008 , Mayer et al. 
|2007| ), the computational expense involved means that these simu- 
lations have, thus far, only barely probed that scales of interest and, 
in doing so, have made restrictive assumptions (typically explicitly 
turning off cooling and/or star formation on small scales); moreover 
they provide only a snapshot at one instant from the parent simula- 
tion - they cannot survey statistical properties of the nuclear region. 
Alternatively some simulations have simply adopted an assumed 
small-scale structure as an initial condition and studied the result- 
ing gas dynamics at small radii (e.g. Schartmann et al. 2009 , Wada 
|& No rman 2002 , Wada et al. 2009). A number of important conclu- 
sions have been drawn from these studies; however, they not only 
bypass the question of the obscuring material origin, but also have 
thus far adopted idealized potentials, without live star formation 
and/or self -gravity of the gas. As such, the appearance and evolu- 
tion of gravitational modes is suppressed. Cuadra et al. ( 2009} and 
|Fukuda et al.| (2000 ) show (albeit in similar idealized studies that 
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neglected star formation and stellar feedback) that when included, 
gravitational torques from self-gravity are an order-of-magnitude 
stronger than hydrodynamic torques from pressure forces or vis- 
cosity; the same conclusions have been reached for intermediate 
(> lOOpc-scale) bars in a large number of hydrodynamic simula- 
( Noguchi|1987|[T988l|Hernquist|1989[|Barnes & Hernquist] 



tions 



1991 



1996[ |Hopkins et al.||2009a] ), and follow from analytic ar- 



guments (see references above and Rice et al. 2005 ; Hopkins & 
|Quataert|20lTa] >. 

Recently, to understand the angular momentum transport re- 
quired for massive BH growth, we have carried out a series of nu- 
merical simulation s of inflow from galactic to BH scales ( Hopkins 
& Quataert^OlOa i. 1 By re-simulating the central regions of galax- 
ies, gas flows can be followed from galactic scales of ~ lOOkpc to 
much smaller radii, with an ultimate spatial resolution < 0. 1 pc. For 
sufficiently gas-rich disky systems, gas inflow continues all the way 
to < 0.1 pc. Near the radius of influence of the BH, the systems be- 
come unstable to the formation of lopsided, eccentric gas+stellar 
disks. This eccentric pattern is the dominant mechanism of an- 
gular momentum transport at < 10 pc, and can lead to accretion 
rates as high as ~ IOMq yr _1 , sufficient to fuel the most luminous 
quasars. In addition, through this process, some of the gas continu- 
ously turns into stars and builds up a nuclear stellar disk. Relics of 
these stellar disks may be evident around BHs in nearby galaxies 
^Hopkins & Quataert|2010b|), such as M31 and NGC4486b dLauer| 
et al.|19931|Tremaine|1995llBender et al.|2005l|Lauer et al.|l996f 



2005 Houghton et al. 2006, Thatte et al. 2000, Debattista et al. 



2006). In this paper, we examine the possibility that the disk that 
drives accretion and accounts for these stellar relics, in its gas-rich 
phase, may in fact be the canonical torus-like obscuration region 
near AGN. If correct, this implies both an a priori understanding 
of torus formation and structure, and an entirely new paradigm in 
which to view the nature of AGN obscuration. 

Specifically, we here perform a first comparison of these hy- 
drodynamic simulations with the observed properties of AGN ob- 
scuration. We focus on dynamical properties and quantities such as 
the column density distribution that can be robustly predicted with- 
out reference to higher-order radiative transfer effects (which will 
be investigated in future work). In §[2] we summarize the properties 
of the numerical simulations, and in § [3] show how they naturally 
form torus-like obscuring structures. In § [4] we consider the scale 
heights and vertical structure of these torii, and examine how this 
can arise independent of stellar feedback from various gravitational 
processes. In § [5] we compare a number of observable dynami- 
cal properties of the predicted torii to nuclear-scale observations of 
AGN. We then in §|6]consider the full column density distribution, 
and in particular how it depends on sub-grid assumptions about the 
dumpiness of the ISM phase structure on un-resolved scales. We 
use this in § [7] to consider the predicted obscured fractions as a 
function of AGN and galaxy properties. Finally, we summarize our 
conclusions and discuss observational tests and future work in §[8] 



2 THE SIMULATIONS 

The simulations described here are from a survey of multi-scale 
"zoom-in" runs which model gas inflows and star formation from 
large galactic scales to sub-pc scales, and have been discussed in a 



1 Movies of these simulations are available at http : //www . cf a . | 
[harvard. edu/ ~phopkins/Site/Movies_zoom. html| 



series of p apers ([Hopkins & Quataert|2010a| |2011a| |2010b[ [Hop] 
|kins|2 010, Hopkins & Quataert 201 1b]). A detailed des cription and 
list of simulations is presented in Hopkins & Quataert (2010a); we 
briefly summarize the salient properties here. 

The simulations were performed with the TreeSPH code 
GADGET-3 (Springel 2005 ); the detailed numerical methods are de- 
scribed there and in Springel & Hernquist (2002); Springel et al. 
(2005b]). The simulations include collisionless stellar disks and 
bulges, dark matter halos, gas, and BHs. For this study, we are in- 
terested in isolating the physics of gas inflow. As a result, we do not 
include explicit models for BH accretion feedback - the BH's only 
dynamical role is via its gravitational influence on scales < lOpc. 

Because of the large dynamic range in both space and time 
needed for the self-consistent simulation of galactic inflows and nu- 
clear disk formation, we use a "zoom-in" re-simulation approach. 
This begins with a large suite of simulations of galaxy-galaxy 
mergers, and isolated bar- (un) stable disks. These simulations have 
0.5 x 10 6 particles, corresponding to a spatial resolution of 50pc. 
These simulations have been described in a series of previous pa- 
pers (|Di Matteo et al.|2005[|Robertson et al.|2006c|b|a[[Cox et al.| 
|2006||Younger et al.|2008[|Hopkins et al.|2009a| >. From this library 
of simulations, we select representative simulations of gas-rich ma- 
jor mergers of Milky-Way mass galaxies (baryonic mass 10 11 Mq), 
and their isolated but bar-unstable analogues, to provide the basis 
for our re- simulations. The dynamics on smaller scales does not de- 
pend critically on the details of the larger-scale dynamics. Rather, 
the small-scale dynamics depends primarily on global parameters 
of the system, such as the total gas mass channeled to the center 
relative to the pre-existing bulge mass. 

Following gas down to the BH accretion disk requires much 
higher spatial resolution than is present in the galaxy-scale simula- 
tions. We begin by selecting snapshots from the galaxy- scale sim- 
ulations at key epochs. In each, we isolate the central 0.1 — 1 kpc 
region, which contains most of the gas that has been driven in from 
large scales. Typically this is about 10 10 Mq of gas, concentrated in 
a roughly exponential profile with a scale length of ~ 0.3 — 0.5 kpc. 
From this mass distribution, we then re-populate the gas in the cen- 
tral regions at much higher resolution, and simulate the dynam- 
ics for several local dynamical times. These simulations involve 
10 6 particles, with a resolution of a few pc and particle masses of 
« 1O 4 M . We have run rsj 50 such re- simulations, corresponding 
to variations in the global system properties, the model of star for- 
mation and feedback, and the exact time in the larger-scale dynam- 
ics at which the re-simulation occurs. Hopkins & Quataert (2010a ) 
present a number of tests of this re-simulation approach and show 
that it is reasonably robust for this problem. This is largely be- 
cause, for gas-rich disky systems, the central ~ 300 pc becomes 
strongly self-gravitating, generating instabilities that dominate the 
subsequent dynamics. 

These initial re- simulations capture the dynamics down to 
~ 10 pc, still insufficient to quantitatively describe accretion onto 
a central BH. We thus repeat our re-simulation process once more, 
using the central ~ 10 — 30 pc of the first re- simulations to initialize 
a new set of even smaller-scale simulations. These typically have 
~ 10 6 — 10 7 particles, a spatial resolution of 0.1 pc, and a particle 
mass « 100 Mq. We carried out ~ 50 such simulations to test the 
robustness of our conclusions and survey the parameter space of 
galaxy properties. These final re- simulations are evolved for ~ 10 7 
years - many dynamical times at 0. 1 pc, but short relative to the dy- 
namical times of the larger- scale parent simulations. We also car- 
ried out a few extremely high-resolution intermediate- scale simula- 
tions, which include ^5x 10 7 particles and resolve structure from 
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~ kpc to rsj 0.3 pc - these are slightly less high-resolution than the 
net effect of our two zoom-ins, but they obviate the need for a sec- 
ond zoom-in iteration and "bridge" the scales of the above simu- 
lation suites. The conclusions from these higher resolution simula- 
tions are identical. 

Our simulations include gas cooling and star formation, with 
gas forming stars at a rate motivated by the observed Ken nicutt] 
( 1998 ) relation. Specifically, we use a star formation rate per unit 
volume p* oc p 3//2 with the normalization chosen so that a Milky- 
way like galaxy has a total star formation rate of about 1 Mq yr -1 . 
Varying the exact slope or normalization of this relation has no 
qualitative effect on our conclusions. However, we caution that 
since we do not resolve the scales of individual bound star-forming 
cores in these simulations, the star formation is probably more uni- 
form over the small radii than it would be in a more realistic ISM 
model. This is unlikely to be important for global properties here, 
but may have important consequences for e.g. detailed radiative 
transfer effects. 

Because we cannot resolve the detailed processes of super- 
novae explosions, stellar winds, and radiative feedback, the effect 
of feedback from stars is crudely modeled with an effective equa- 
tion of state (Springel & Hernquist 2003}. In this approach, feed- 
back is assumed to generate a non-thermal (turbulent, in reality) 
sound speed that depends on the local star formation rate, and thus 
the gas density. Hopkins & Quataert (2010a ) describe in detail the 
effects of different subgrid ISM sound speeds on angular momen- 
tum transport and inflow rates, and argue that observations favor ef- 
fective turbulent speeds of ~ 10 — 50 kms -1 for densities ~ 1 — 10 5 
cm -3 , respectively. But because the real physics and their effects 
are uncertain, it is important to vary this prescription and deter- 
mine which of our conclusions are sensitive to the assumed sub- 
grid properties. 

Within the context of this model, we can interpolate between 
two extremes using a parameter q eos . At one end, the gas has an ef- 
fective sound speed of lOkms -1 , motivated by, e.g., the observed 
turbulent velocity in atomic gas in nearby spirals or the sound speed 
of low density photo-ionized gas; this is the "no-feedback" case 
with geos = 0. 2 This is broadly similar to what is assumed in 



Bour- 



|naud et ah] ( |2007| >; |Teyssier et al.| ( |2010| ). The opposite extreme, 
^eos = 1, represents the "maximal feedback" model of Springel 
|et al.|p005b) ; in this case, 100% of the energy from supernovae 
is assumed to stir up the ISM. This equation of state is substantially 
stiff er, with effective sound speeds as high as ~ 200 kms -1 . This is 
qualitatively similar to the near-adiabatic equations-of- state in the 
BH accretion studies of |Mayer et~aT] ( |2007| >; |Dotti et al.| (2009) . 
The sound speed at scales we consider cannot meaningfully be 
much larger than this, since it is similar to the circular/escape veloc- 
ity. By varying q eos , we examine a spectrum of intermediate cases: 
for example, equations of state similar to the "starburst" model in 
Klessen et al. ( 2007 ) or the sub-GMC equation of state in |Spaans] 



& Silk| ( |2005| >. Most of our suite of simulations focuses on a wide 
range of sub-grid sound speeds ^ 20 — 100 kms -1 , motivated by a 
variety of observations of dense, star forming regions both locally 
and at high redshift (Dowries & Solomon 1998, Bryant & Scoville 



2 This is still a non-trivial dispersion at large radii in galaxy disks. At the 
scales we focus on here, however, this corresponds to sounds speeds far 
below the circular velocity, and Jeans masses ~ 100 Mq, our resolution 
limit. As such, allowing cooling to even lower temperatures = 10 K makes 
no difference beyond the g e os = case. 



1999; |Forster Schreiber eFaTl[2006l |Tono et al.|20fi7) , and recent 
numerical simulations ( Hopkins et al. 2011a). 

Within this range, we found little difference in the physics of 
angular momentum transport or in the resulting accretion rates, gas 
masses, etc. on the scales we consider ( Hopkins & Quataert 2 010a| ). 
More detailed comparison with the explicit stellar feedback models 
presented in Hopkins et al. (20 11a|c|b| ) will be the subject of future 
work. Here, we will focus on the effects on the obscuring gas near 
the BH. Because we are not explicitly accounting for or resolving 
feedback processes, we do not expect these models to accurately 
reflect the detailed dynamics of gas in response to strong feedback. 
Rather, we wish to use our suite of simulations to identify behavior 
that is robust to the effective pressure or turbulent sound speed of 
the gas - i.e. to identify robust aspects of the system that are present 
even without feedback such as stellar winds. 



3 FORMATION OF THE TORUS 



Hopkins & Quataert (2010a) show that when large-scale inflows 
are sufficient, the buildup of gas in the central regions of the galaxy 
triggers a cascade of secondary instabilities, that drive rapid inflows 
to still smaller radii and ultimately onto the BH. Around the BH 
radius of influence, these instabilities generically take the form of 
an m — 1 mode - a thick, eccentric, slowly precessing gas+stellar 
disk, in which the eccentric stellar pattern torques strongly on the 
gas, inducing shocks and inflows. The disk can then propagate gas 
inflows and the m — 1 pattern down to small radii < 0.1 pc, where 
it transitions to a traditional alpha-disk. This should be generic to 
any quasi-Keplerian potential in a dissipative system with shocks 
( [Hopkins &^uataert|20fTal|Hopkinsp0T0] ). 

Figure [T] shows some illustrative examples of the nuclear gas 
disks that form around the BH radius of influence in our simula- 
tions. We plot gas surface density maps, with color encoding the 
gas effective sound speed, from scales of > lOkpc to < lpc. The 
initial large-scale simulation in this case is a fairly gas-rich ma- 
jor merger of two ~ L* galaxies (with initial bulges of mass 1/3 
the disk mass and BHs of mass 10 7 Mq). The zoom-in simulations 
were carried out just after the coalescence of the two nuclei, which 
is near the peak of star formation activity, but when the system is 
still quite gas rich. 3 

We both show the global structure, and edge-on (R, z) disk. 
The scales shown include the BH radius of influence, about lOpc 
in these galaxies. In the face-on projection, the m — 1 modes that 
form at these scales are clearly evident. They drive large torques on 
the gas, driving inflow into <C 0. 1 pc at accretion rates as high as 
IOMq yr -1 in these simulations, sufficient to power the most lumi- 
nous quasars (see Figures 5 & 13 in Hopkins & Quataert (2010a)). 

Here, however, we note the broad resemblance of these nu- 
clear disks to the canonical AGN "torus." The disks are thick, with 
characteristic scale ^ 0.1 — lOpc, gas masses ~ M B h, and scale 
heights of order unity. Of course, unlike in toy models of the torus, 



3 The specific properties of each simulation are given in |Hop-| 
|kins & Quataert] {20 10a), those shown here are (top-to-bottom, 
left-to-right): Nf8hlc0thin, NfShlclthin, Nf8hlclqs, Nf8hlcldens, 
Nf8hlc0 (left); Nf8hlclICs, Nf3hlclmid, Nf2h2b2, Nf8h2b2, 
Nf8h2b4 (right). They have (respectively) initial gas fractions 
/ g as ~ 0.5,0.6,0.8,0.8,0.8,0.75,0.26,0.20,0.8,0.8; BH mass ~ 



3 x 10 7 Mq and disk mass ' 
10 7 Mq inside lOpc, 
35, 20, 40, 50, 10, 40, 30, 25, 25, 20 kms" 



< 1.2,1.7,3.0,8.1,0.25,1.7,4.6,7.0,3.5,0.5 x 
and sub-grid sound speeds c s ~ 
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Figure 1. The face-on (x-y) and edge-on (x-z) disk structure of the nuclear disks in several representative simulations. Scale is the same in all panels (lower 
right). Each example is a simulation of the central ~ lOOpc of galaxy nuclei with different initial large-scale galaxy properties (inflow-to-BH mass ratios, 
gas fractions, and treatments of stellar feedback; details in text). Intensity encodes gas surface density (increasing from Nh < 10 21 cm -2 to Nh ^ 10 25 cm -2 . 
Colors encode the absolute star formation rate of the gas (increasing from blue to red/yellow). The formation of a lopsided, gas-rich disk is ubiquitous. Regions 
where gas shocks (edges in this image) dissipate energy, leading to rapid gas inflow. Viewed edge-on, the disks are all thick, with columns > 10 22 cm -2 to 
h/R ~ unity. 



the gas is part of a continuous distribution at all radii, and its struc- 
ture is non-trivial. 



4 VERTICAL STRUCTURE: DEPENDENCE ON 
STELLAR FEEDBACK 

4.1 Overview 

The major input parameter of our models is the parameterization 
of the effects of stellar feedback on the ISM. This is accomplished, 
here, with the parameter q CO s described in § [2] that allows us to 
interpolate between a feedback-free ISM and one with large non- 
thermal internal gas velocities and pressures driven by stellar feed- 
back. 

The so-called torus is defined largely by its vertical struc- 
ture, which determines the obscured fractions. To the extent that 
the amount of turbulent velocity and pressure support in the simu- 
lated gas is defined by a sub-resolution model, we must ask whether 
the vertical structure we see in our simulations is entirely a conse- 
quence of our model inputs, or whether there are robust statements 
and predictions we can make. 

We therefore consider the vertical structure in detail, in a spe- 
cific su rvey of q eos . This survey (Nf8h2b4q in Hopkins & Quataert 
2010a ) is a typical, canonical set of conditions (3 x 10 7 Mq BH, 



with disk-to-BH mass ratio of a few inside ~ lOOpc initially, and 
initial gas fraction ~ 50%, typical of the simulations in Figure [TJ. 
We re-simulate the identical cases, but with q eos = 0.0, 0.018, 0.06, 
0.10, 0.12, 0.15, 0.21, 0.25, 0.35, 0.60, 1.0. The spacing in q eos 
is chosen such that the implied turbulent gas sound speeds are 
spaced over roughly equal logarithmic intervals from the minimum 
<7eos = floor (lOkrns -1 ) to the maximum q eos = 1 value (which is 
density dependent, but ~ lOOkrns -1 at range of interest). 

Figure [2] shows the edge-on (R, z) gas structure, as a function 
of ^eos- A few generic features stand out. The disks are generally 
thick. At the smallest radii (< 0.1 — lpc), they eventually become 
thin, since the gravity from the BH becomes arbitrarily strong. This 
gives a torus-like morphology. Flares (discussed below) and lopsid- 
edness (reflecting the lopsided disk mode on these scales) are not 
uncommon. As a function of q e0 s, we see unsurprisingly that the 
gas distribution becomes more smooth and vertically extended at 
higher q eos . For q CO s > 0.4, the system is no longer really a vertically 
supported disk, but spherical - however, as discussed in Hopkins 
|& Qu ataert (20 10a| ), this is likely an unrealistically large implicit 
feedback efficiency. 

The most surprising thing about Figure [2] is how little change 
there is as a function of For q eos ~ — 0.35, there is a factor of 
~ 5 change in c s , which leads to a naive expectation of a factor of 
^5 — 25 change in h/R. We see much weaker variation. We now 
consider this quantitatively. 
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Figure 2. Edge-on gas mass distribution plotted in cylindrical coordinates (R, z) to highlight the torus-like structure of the disk, with intensity and color as 
Figure^ The simulations shown here are a survey of q CO s, which determines the effective (sub-grid) pressure support of the ISM: the (mass-averaged) effective 
sub-grid sound speed c s is labeled in each panel. Figure^ for some of our survey of q eos . As expected, the systems become more puffed-up with increasing 
g e os (sub-grid c s ), and for q eos > 0.4 they are nearly spherical. But at small q eos , the scale heights do not decrease as rapidly as oc c s , but approach some 
asymptotic minimum. 



4.2 General Expectations 

To inform our comparisons, consider a simple smooth, isothermal 
system, in which the self-gravity of the gas is negligible (i.e. the 
potential is dominated by the BH, stars, and/or dark matter). The 
equation of vertical hydrostatic equilibrium 

OP 2 1 dp d$ 
oz p oz oz 
then has the trivial solution 

p(R,z) = P o(R)exp{c7 2 [$(R,0)-<S>(R,z)]} . (2) 

For large z/R, this depends on the specific form of <£, and so on 
the details of the global mass distribution. However, if the disk is 
thin, i.e. most of the mass is at z/R <C 1, then this has a particularly 
simple expression. For any background spherical mass distribution, 
we have d<S>/dz (d<S>/dr) (dr/dz) (V 2 /r) (z/r), where r 2 
R 2 + z 2 and V 2 = GM enc ( < r) jr. So for z 0) - $ (fl, z) « 

GMenc(< R)R~ 2 z/2 « V 2 (z/R) 2 /2. 

Together, this gives the especially simple solution for the den- 
sity for a quasi- spherical potential: 

pOM«p (fl)expj-i (£j | (3) 

h l= cs_ ( GM enc (<R) Y 1/2 

R ~ V c ~ ° s \ R ) (4) 

Of course, the c s here does not need to be thermal. Non- 
thermal pressure sources such as turbulent motions will have the 
same effect. So for comparison with simulations we should take 



c s — > c Z) es, where c\ eff = c 2 + o\ includes both the thermal and/or 
sub-resolution effective sound speed (c s ) and resolved turbulent 
vertical motions (<j z ). 

Figure [5] compares this expectation for p(z) as a function of 
£z,eff/Vc to the actual vertical mass distribution measured in nar- 
row radial annuli from ~ 1 — lOpc. We use the full c z , e ff as defined 
above. The distributions are reasonably described by the above 
scalings, A gaussian core is typical, with a slightly broader (often 
more exponential) distribution at high-z. Remember that at suffi- 
ciently large \z\, the correct solution involves the full potential; if 
we account for this more accurately, we see similar agreement. The 
important point is that the gas does appear to be in vertical equilib- 
rium. 



4.3 Gravitational Support 

Given the gas dispersion c z?e ff in the simulations, the vertical struc- 
ture is what we would expect. But are these dispersions primarily 
sub-resolution (set by the model), thermal, or gravitational in ori- 
gin? 

Figure [4] again considers the vertical gas profile, but in our 
survey of different q e0 s- We know from Figure[3]that accounting for 
the full gas motion explains the observed scale heights. Therefore 
here we compare the expectation if the gas motions were purely 
thermal and/or sub-resolution - i.e. c ze ff = c s , where c s is the sound 
speed and is dominated by the sub-resolution turbulent effective c s 
(since the explicit cooling time of the gas is ~ 10 4 times shorter 
than its dynamical time). 

In the higher-^eos (higher effective c s ) simulations, this ex- 
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Figure 3. Vertical gas mass distribution at different radii. Panels show the integrated gas mass fraction (M gas (> \z\)/M gas ) above a given height |z|, for gas in 
a narrow radial annulus about R (each R labeled). The heights are normalized to the expected scale height zo = c e ff/K, where c e ff = (c 2 s + a\ ) 1//2 , where c s 
is the sub-grid velocity dispersion (plus thermal sound speed) and cr z is the resolved gas velocity dispersion. Each line is a different simulation (with varied 
initial gas fraction, disk and BH mass, sub-grid equation of state, and star formation laws); shown at a randomly-chosen time near the peak of the inflow onto 
the BH (but behavior is similar over the entire duration of the simulations). Thick dashed black line is the simple Gaussian expectation for an isothermal gas 
disk with weak self-gravity in vertical equilibrium (Equation[3j. Bottom Right: The scale height h (best-fit dispersion zo, fitting the vertical gas distribution at 
each radius to a Gaussian), as a function of radius, for each simulation. The scale heights are significant, and the vertical behavior approximately follows the 
linear expectation at these radii, if the full vertical dispersions are included. 
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R [pc] 



R [pc] 



R [pc] 



R [pc] 



Figure 4. Top: Vertical gas distribution, in our survey of q e0 s- Each panel shows a different simulation in the q QOS survey. Each line shows the vertical gas 
distribution at a different radius (as labeled). Here, the x-axis is scaled only by c s /V c , the scale height expected if the sub-resolution (feedback-driven) velocity 
dispersions were dominant (as compared to c e ff). Bottom: Gas scale height h versus radius (solid black lines). Dotted blue lines show the expected height if 
just the subgrid velocities were present. Dashed red lines compare the expected height including the full resolved dispersion c e ff. In large-g e os systems, the 
implicit feedback dominates the vertical support. But the scale heights in low-g e os systems are supported by large resolved turbulent vertical velocities, despite 
the lack of feedback. There is some non-feedback dispersion source in these systems. 
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plains most of the pressure support, i.e. the resolved turbulent dis- 
persion a z <C c s (sub — grid). But at low q e0 s, the scale heights and 
c z ,e& do not drop nearly as quickly as the sub-grid c s alone. There 
is some non-thermal, resolved gravitational process giving rise to 
minimum scale heights and vertical dispersions. 

What, then, dominates the effective vertical "heating" in the 
torus region? 

4. 3. 1 Clump-Clump Encounters 

It has been proposed that two-body scattering between dense 
molecular clumps in the gas could maintain the observed scale 
heights jKrolik & Begelman|1988[[Nayakshin & King|2007[|Hobbs| 
|et al.||2010| ). However, we find these effects are negligible in our 
simulations. 

Consider clumps within the plane of a disk. Scattering a clump 
to large V z ~ V c requires both (a) an encounter between two clumps 
with relative velocity > V c , and (b) an encounter within an impact 
parameter b such that GM c \/bV c > V c . The mean time per clump 
between such encounters is just r _1 ~ f(V c /a)n c \b 2 V c , where 
n c \ is the volume density of clumps and f(V c /a) is the fraction 
of the clumps moving on orbits with large non-circular motions 
(| V — V c \ > V c ). If the system is sufficiently thin such that b > h, the 
disk thickness, then this becomes r _1 ~ f(V c /a)dN c \/dAbhV c . 
Using n c \ = pg as /M c \ — Y, gas /hM C h the maximum b above, and 
Vc rsj GM enc /r, this can be written 



f(VcM 



Menc 
M 2 as 



(1 + gtfcl) 



(5) 



where Q is the Toomre Q ~ (h/R) (M ga% /M enc )~ l and N c \ the total 
number of clumps, and the expression shown interpolates between 
the extremely-thin and thick-disk cases. For a Maxwellian veloc- 
ity distribution, f{V c /cr) ~ exp{-(V c /cr) 2 /2}. Since both/(V c /cr) 
and Mg as /Menc ~ Mgas/MsH are small at this radius, collisions re- 
quire many dynamical times. But any induced vertical heating will 
relax away in just a single or couple dynamical times, since the 
cooling time is much shorter than the dynamical time. So without 
continuous energy input to drive large dispersions - which is es- 
sentially the problem we wished to solve in the first place - this 
mechanism fails. 

Moreover, if star formation occurs with some efficiency rela- 
tive to the dynamical time (p* ~ ep y/Gp, with e ~ 1 — 10%), then 
using the fact that any clump must have p > M enc /R 3 to avoid tidal 
destruction, clump-clump gas heating must occur faster than the 
gas exhaustion timescale in a clump, requiring 



(6) 



Even for a ~ V c (which begs the question) and an extremely thin 
disk QN C \ < 1, this requires M gas /M enc ^> 0. 1, which is not satisfied 
at the inner radii < lOpc. 

4.3.2 Twists and Misalignment 

Another possibility is that large covering factors are maintained 
by virtue of the fact that the nuclear disk is mis-aligned with the 
larger- scale inflow/bar/disk. This is particularly interesting because 
observations find relatively little correlation between the axes of 
AGN (traced by jets or the torus) and the inclination of the host 
galax y (e.g. |Keel||19801 |Lawrence & Elvis|[T982l |Schmitt et al. 
17997) |Simcoe et al.||1997| |Kinney et al.||2000[ |Gallimore et al' 
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Figure 5. Contribution of disk mis-alignments and twists (as opposed to 
disk thickness) to obscuration. For each simulation in the sample shown 
in Figure [3] we calculate the (time-averaged) fraction of sightlines towards 
the BH which would be obscured (Nh > 10 22 cm -2 ) assuming the disk is 
razor- thin at every radial annulus Ar (using the angular momentum vector 
of gas in A r to define the disk plane). If the disk plane (angular momentum 
direction) were constant with radius, this would be zero, but if the plane 
tilts as a function of radius, it can be non-zero (a 180° flip at some radius, 
maintained for the duration of the simulation, would give an obscured frac- 
tion of unity). Although there are significant mis-alignments (see [Hopkins 
|et al.|20lTd) , and in a few systems they can account for obscured fractions 
> 50%, they would only give an integrated obscured fraction of ~ 25% if 
the disk were thin at all radii. And even simulated cases with no "twists" 
still have large obscured fractions and h/R in Figure[3] The torus must ac- 
tually be thick to match observations - so some process must explain large 
scale-heights in the gas. 



2006 ;|Zhang et al.|[20Q9| but see also Mai olino & R ieke 1995, 
Shen et al.||2010| and references therein). In a companion paper, 
Hop kins et al.| ( 20TTd| ), we show that this lack of alignment is re- 
produced in our simulations owing to two processes. First, occa- 
sionally the central gas supply is strongly influenced by a single 
or couple large clumps that form at large radii, fragment and sink, 
realigning the central angular momentum vector (see e.g. |Nayak-| 
shin & King 2007 ). Examples of this have also been seen in cos- 
mological zoom-in simulations ( [Levine et al.|2010| ). Second, even 
in smooth flows, the secondary and tertiary gravitational instabili- 
ties will tend to de-couple their angular momentum from the pri- 
mary (external) bar/spiral structure, and semi-chaotically precess 
or tumble in three dimensions (see Heller et al. 2001 ; Shlos man & 
Helle r|2002(|El-Zant & Shlos man|2003{ Maciejewski & Athanas- 
|soula|2008||Englmaier & Shlos man|2004) . 

In Figure|5]we show how this can contribute to obscuration. It 
is straightforward to measure the axis of angular momentum of the 
disk in a radial annulus, and define the corresponding inclination 
G(r) (relative to the initial, uniform angular momentum axis of the 
entire initial disk). We also know the mass AM gas (r) enclosed in 
each annulus; we can simply integrate along all sightlines towards 
the BH, assuming the mass is in an axisymmetric razor-thin disk 
with inclination O(r), to obtain the column density £ gas along each 
sightline. If a sightline is covered by the disk at some radius, it is 
"obscured." 4 



4 Technically, we require a column that translates to Nh > 10 22 cm -2 , but 
because of our razor-thin assumption, this is almost identical to being cov- 
ered by the disk. 
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We consider these assumptions because they effectively de- 
fine a minimum obscured fraction stemming purely from twists 
and misalignments. This fraction can be considerable, but there is 
a broad range in different simulations - many systems have only 
< 20% covering fractions, but there is a long tail towards near full 
covering (anti- alignment of the central and outer disks). Integrated 
over all simulations and snapshots, the average covering fraction is 
w 25%. A warped or twisted disk can therefore yield large covering 
angles towards the BH even when the disk itself is thin. 

However, this is not the full story. First of all, the covering 
fractions of « 25% are still significantly lower than the total cover- 
ing fraction of obscuration in the simulations, by a factor of at least 
~ 2. Moreover the cases with weak twists (<^C 20% covering in Fig- 
ure^ still exhibit large obscured fractions and thick disks. The key 
point is that the vertical density distribution in Figures |3|4| shows 
that we must explain the actual thickness, not just the orientation of 
the disks. This is true for observations as well - empirical model- 
ing of the hot dust continua indicates that the obscuring region must 
be geometrically thick, not just a misaligned larger-scale thin disk 
(e.g. |Deo et aL]|2009| and references therein). A time-dependent 
twist can, in principle, "pump up" vertical motions, but fast cooling 
times make it difficult to sustain a large scale height anywhere ex- 
cept close to the location of the twist (where the pumping occurs). 
Some mechanism that pumps vertical motion throughout the disk, 
on a timescale comparable to the local dynamical time, is required. 



4.3.3 Bending Modes 

Bending modes can provide an efficient channel for "heating" the 
torus. Their behavior is particularly interesting in response to "slow 
modes" in a quasi- Keplerian potential. Consider a general bending 
mode 



h(R, 0, t) - H(R) exp {i (uj b t - m b 0)} 
H{r) = ho(R) exp ji J k b {R')&R' 



(V) 
(8) 



in a system that includes some quasi- spherical component 
(BH+bulge+halo) and a thin disk with surface density E</, angular 
(vertical) frequency Q (is), and velocity dispersions in the radial, 
azimuthal, and vertical directions a r , a^, a z . The value fa is the 
radial wavenumber of the bending mode, and m b is its azimuthal 
wavenumber. In the WKB regime, if a 2 > a 2 ^, the dispersion rela- 
tion can be written 



(cj b -m b Q) 2 = z/ + 27rG£j \k b \ + (a 2 - a 2 ) k\ 



(9) 



flKulsrud & MarklfTWOj |Kulsrud et aLlfT97T] pa^|[T97T| |PolTl 
achenko|!977| >. 5 

If 07 is sufficiently large, the system is vulnerable to the so- 
called "firehose" instability and bending modes will be self-excited. 
However, it is unusual to see such large a r (and a 2 > a^) in disks. 
Even with large a r , the fact that a r < V c (r) for any meaningful 
"disk" means that usually, when the self-gravity of the disk is small 
compared to the background potential, the system is stable. And in 
even in self-gravitating disks with large a r , it typically takes only a 
small v z to stabilize them, so the induced h/R is not large. 

5 Note that the a 2 that appears in Equation's not technically a dispersion 
(that being defined (v 2 ) — (v) 2 ), but the mean (v 2 ). Thus streaming/bulk 
motion in the radial direction is affected just as much as random motions 
about some mean v r (important for our purposes, since gas parcels being 
collisional tend to move in coherent streaming motion). 



However, consider the special case of interest here, where 
the disk is quasi- Keplerian and has a large lopsided mode driv- 
ing accretion. The system is (initially) a thin disk in the quasi- 
Keplerian potential of a BH - i.e. to lowest order, the parame- 
ters are those of a pure Keplerian potential, with some correction 
terms that scale with Mj/Mbh of O(e) < 1. As discussed above, 
and in previous works, the disk develops a gravitational instabil- 
ity in the disk plane (the standard density waves of spiral/bar/etc. 
modes), which we can describe by e.g. the perturbed density field 
Ei(fl,f) = |a|E (fl) ^{i(J R k p (R')dR' ^uj p t-m p 0)}.Here \a\ 
is the effective mode amplitude in the density field at R, and the 
properties uj p ,m p , and k p refer to the frequencies and wavenumbers 
of this, in-plane mode (independent from the uo b ,m b , and fa of the 
bending mode). The fact that the potential is quasi-Keplerian, i.e. 
has « k, favors (and supports for long periods of time) global, 
"slow" m = 1 modes - modes with m p — \, lu p ~ eQ <^ Q, and 
\k p R\ ~ 1. These are the lopsided/eccentric modes that we see 
above. The potential of the BH+disk system is 



GM Bn 



and it is useful to define the parameter 



2Q 



1 

~2Q 



2d_ + d^_ 
r dr dr 2 



(10) 



(11) 



To first order in e, then, the WKB dispersion relation of such 
modes in a cold (c s <C V c ) disk is 



uj p — zu -\-itGYid \k p \ Q 



(12) 



(Tremaine 2001). The equations of motion for the perturbed veloc- 
ity v = R Q cj) + v r R + <f> become, at this order, 



V r - 



2(uj p - 



d&i 2$i 
dr r 



-n\k P r (i3) 



V0 = ^ V r 

where we have used the WKB relation <3>i 
Since a 2 = (|v r | 2 ) and V c = QR, this becomes just 



(14) 



-27rG|^| _1 Ei. 



a r =\a\\k p R\- 1 Vc 



(15) 

with a 2 =4(7^. And recall for our simulations, the magnitude \a\ 
in the disk is order unity during the active phases of BH growth 
( [Hopkins & Quataert|2010a| >. 

This is the important point - because of cancellations that oc- 
cur (essentially the entire system is near-resonance), the radial in- 
duced velocities from the m — 1 mode are quite large, ~ V c , inde- 
pendent of how small the ratio Mj /Mbh may be. 

Now return to the dispersion relation for bending modes. For 
the non-disk (Keplerian) part of the potential, v 2 — Q 2 . To lead- 
ing (zeroth) order in e ^ Mj/Mbh, then, the dispersion relation be- 
comes 

2 



1 + 



\a\ 2 \k p R\ 



Mr 



(16) 



where we have defined h/R — cj z /V c . 

Recall, lall&p/?! -1 ^ 1 for the lopsided disk mode, and 
\faR\ ^> 1 in the WKB limit - thus, whenever the disk is thin, the 
radial motions induced by the m p — 1 eccentric disk excite bend- 
ing modes. These bending modes will grow on the local dynamical 
timescale, since uo b ~ f2. Compare the slow modes in the plane, 
which have u p ~ e Q and can be long-lived. 

The bending mode, of course, creates some non-trivial vertical 
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Figure 6. Comparison of the scale-height h (as in Figure [3} to the m = 1 
eccentric disk mode amplitude \a\ = at different radii (each point 

samples radii evenly in logR from R = 1 — lOpc), for different simula- 
tions (each color denotes a different simulation, near its peak activity). All 
simulations shown are "cold" (i.e. have weak stellar feedback assumed in 
their sub-resolution prescription, specifically q e0 s < 0.1), so that the scale 
height is dominated by resolved turbulent vertical motions, not the feedback 
model input. In these cases, there is a correlation of the form h/R ~ \a\, 
corresponding to the prediction from an h /R pumped-up by bending modes, 
themselves excited by radial motions from the in-plane eccentric mode. The 
torus height can be continuously sustained by exchange of energy from the 
eccentric/lopsided disk mode that powers the accretion onto the BH. 



motion. The growth of the mode will saturate when it drives a a z 
sufficiently large so as to reach the marginal stability condition of 
Equation 16 (Im(u;&) = 0), which we can write as 



•Ml" 1 (W\ 2 \k P R\- 2 \hR\ 2 -l) (17) 

(18) 



where the second equality uses \foR\ > 1. In short, tightly- wound 
bending modes arise, and saturate v z at a large fraction of V c , such 
that h/R is driven to order unity wherever the eccentric mode per- 
sists, independent of the degree of self-gravity of the disk! 

In Figure [6] we check whether this prediction at all describes 
our simulations. We compare the scale height h/R to the measured 
mode amplitude \a\ of the in-plane m = 1 mode, at a random time 
during the active phase, for each of a subset of our simulations. We 
chose only the simulations for which q eos < 0.1, where we can con- 
firm that the sub-grid assumed c s does not dominate c e & or the ver- 
tical scale height on the scales we measure (see Figure[4]). We sam- 
ple both quantities at even intervals in \ogR from R = 0.3 — lOpc. 
There is, unsurprisingly, large scatter, but a correlation is significant 
at > 3 a and consistent with h/R ~ \a\ over most of the simulated 
range. That the relation is not exactly linear at the high-h/R end 
and shows considerable scatter is expected, both because of contri- 
butions from kb and k p in the derivations above, non-linear effects 
(especially at h/R and/or \a\ > 0.1), and some non-zero support 
from c s . But it is quite unlikely that this relation would arise acci- 
dentally - after all, for otherwise equal properties, a lower-/*//? disk 
is actually more gravitationally unstable, so if anything we would 
naively expect the inverse of the observed correlation. 



5 BASIC DYNAMICAL PROPERTIES OF THE "TORUS" 

Thus far, we've focused on the origin of torus structural properties 
in simulations. We now examine these properties in more detail and 
compare to observations. Figure [7] shows a number of (azimuthally 
averaged) properties of the nuclear gas, as a function of radius. We 
plot the gas surface mass density, gas fraction, SFR, vertical gas 
velocity dispersion, and gas inflow rate M (here defined so posi- 
tive is inflow). The velocity dispersion includes both resolved and 
sub-resolution components, i.e. c% s = c 2 s + of , where c s is the sub- 
grid implied sound speed (plus any thermal components) and a z is 
resolved vertical dispersion. 

We show this for our suite of simulations from Figure [4] in 
which we systematically vary the sub-grid equation of state (via 
the parameter q e0 s)- For each, we select a random snapshot near the 
peak of inflow activity. Because the global properties - gas density 
profiles, inflow rates, circular velocities, etc - are primarily set by 
global gravitational torques (see Hopkins & Quataert 2011a), the 
parameter q eos does not appear have a dramatic qualitative effect on 
these properties. The primary effect is to determine the efficiency 
of fragmentation, which in turn changes the variability and global 
efficiency of star formation and gas exhaustion. If we consider the 
wider range of simulations shown in Figure|3] which vary the initial 
gas fractions, bulge-to-disk, and BH-to-disk mass ratios, we find a 
similar range in the predicted properties. 

In more detail, Hopkins & Quataert (201 lb) show that the sur- 
face density profiles that arise are a natural consequence of the 
dynamics of tidal torques from the m = 1 lopsided disk instabili- 
ties. Specifically, the perturbation dynamics set a robust range of 
"quasi-equilibrium" profiles in which the gas mass density remains 
quasi- steady state over the active phase so long as there is suffi- 
cient initial inflow to trigger the process. If the profile is a power 
law S oc R~ v , then this range is 1/2 < 77 < 1, similar to that seen 
in "cuspy" ellipticals. 

The SFR surface density follows simply from the assumed lo- 
cal relation between star formation efficiency and dynamical time: 
in the simulations, p* oc p 3 ^ 2 . Competition between gas inflows and 
SF sets the gas fractions, although these evolve significantly via de- 
pletion. 

There are some observations to which we can compare. Wa- 
ter masers have been observed and used to map the inner disk 
structure around AGN in a few nearby galaxies ( Greenhill et al 



[T997 2003 ; Br aatz et al.|2004||Henkel etall2 005 ; Kondratk o et al 
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2006b a 2008). These are sensitive to densities ~ 10' 
(typically ~ 0.1 — lpc). At larger radii, interferometry has also 
been used to image the molecular and HI gas in the nuclei of 
some nearby systems (Lonsdale et al. 2003 ; Schinnerer et al. 2000 
Combes et al. 2004 ; Garcfa-Burillo et al. 2005 ; Schinnerer 



2008). Complemented with adaptive-optics imaging of nearby nu- 



clei, this gives constraints on the gas+stellar dynamics, and in- 
formation on the star formation history ([Kuntschner et al. 11200 1| 
IDavies et al.|2006| [Sanchez et al.|20061|Davies et al.|2007HHicksl 
|et al.|2009) . 

We compile these observations and compare to our simula- 
tions in Figure [7] Most of the observed systems have BHs with 
broadly similar masses to our ~3x 10 7 Mq. We plot the obser- 
vations at all radii available. The maser observations are shown as 
points with error bars for resolved properties of disks outside the 
minimum radius enclosing the BH. The larger-scale surface den- 
sities mapped from the gas velocity fields with VLBI are shown 
as solid lines. For constraints involving stars (gas fractions, SFR), 
the VLB 1+ AO constraints are shown as diamonds, at the minimum 
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Figure 7. Azimuthally-averaged nuclear disk properties versus radius. Each dotted line is a simulation with different galaxy and inflow properties, but one 
where some nuclear inflow is excited. We chose a random time near the peak of activity for each simulation to show here, but the results are similar over the 
entire active phase of each. Here for a set of simulations with identical initial conditions but varied q e0 s, the parameter describing our effective equation of state 
and stellar feedback model. Lines range from black (g e os = 0; effective c s = 10km s _1 ) to red (g e os = 1; effective c s = lOOkrns -1 ), evenly spaced in logc 5 . 
Top Left: Surface mass density profiles. Points with error bars show constraints from AGN maser disks (NGC 3079; Kondratko et al. 2005 magenta), (NGC 
3393; |Kondratko et al.|2008| dark green), (NGC 1068; |Lodato & Bertin|2003| pink), (Circinus; |Greenhill et al.|2003| light green). Solid lines show constraints 
from adaptive-optics (AO) measurements of A GN (NGC 1068, 10 97, 3227, 3783, 4051, 4 151, 6814, 7469, Circinus; in violet, green, cyan, blue, dark green, 
magenta, red, orange, dark blue, respectively) (Davies et al. 2007) and Hic ks et al.|{2009) . Top Center: Gas fraction (M gas (< R)/(M* (< R) +M gas (< R))). 
Diamonds are the AO systems (same color styles), at the resolution limits for complementary constraints used to derive / gas . Top Right: SFR surface density. 
For each AO system where a measurement is available, two points are shown. The first is the current SFR (small points), second is the peak SFR (larger points), 
both estimated from the fits to the SFR history inside the minimum and maximum observed radii in |Davies et al.|p007) . Bottom Left: Integrated SFR(< R). 
Bottom Center: Vertical velocity dispersion of the gas. Dotted lines show both the resolved and sub-grid assumed dispersions (added in quadrature). Bottom 
Right: Instantaneous inflow rate through R, generated by gravitational torques from the eccentric disk structure. 



resolved radii of the AO observations. The nuclear SF history is 
modeled for several cases in Davies et al. ( 2007); we show their es- 
timated current SFR both at the innermost radii where stellar light 
is measured and at the outer radii where the integrated light is used 
to determine the SFH. We also show their estimated maximum SFR 
of each observed burst from the fitted SFH within the observed ra- 
dius. In all cases the observations broadly bracket the simulations, 
albeit with larger uncertainties in / gas and the SFH. 

Of course, since these properties all scale with the dynami- 
cal properties of the system, they are all mutually correlated. A 
Kennicutt- Schmidt type law similar to that observed (Esf oc S gas ; 
for the nuclear-scale observations see |Hicks et al.|2009] > is effec- 
tively built into our simulations by sub-grid assumption. 6 We have 
discussed extensively the gravitational origin of the dispersion (a). 
But both a and £ are related to V c , for obvious dynamical reasons, 
and increase at smaller radii and/or in more massive/dense systems. 
And E is tied to Esf via the Kennicutt relation. We therefore pre- 
dict a relation between Esf and a for purely gravitational dynamic 

6 Both the observed and simulated Kennicutt-type laws appear to have an 
index closer to 77 ~ 1.7 rather than the canonical 77 ~ 1 .4. In the simulations, 
this is because we assume a local oc p 3 / 2 , and for the simple case of a 
gas disk contracting at constant h/R this predicts 77 = 1.75. 



reasons. In the past, such a correlation has been interpreted as evi- 
dence of stellar feedback driving the observed dispersions - we find 
this may not be necessary. 



6 THE COLUMN DENSITY DISTRIBUTION: TO CLUMP 
OR NOT TO CLUMP? 

Thus far, all of our analysis has concerned global properties of the 
simulated torii, which we have reason to believe should be robust 
to the exact micro- structure of gas on unresolved scales. However, 
sub-resolution structure can be important in calculating the column 
densities observed towards the BH. We therefore consider this now 
with two simple sub-resolution models. 



6.1 The No- Substructure Case: Smooth Torii 

One extreme is trivial: we simply take the gas distribution exactly 
as-is from the simulations, without any assumed sub-grid substruc- 
ture. The column density along a given line-of-sight at each time 
can then be simply determined (following Hopkins et al. 2005c ). 
We generate ~ 1000 radial lines-of-sight (rays) uniformly spaced 
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in solid angle and with its origin at the BH, and integrate the line- 
of- sight density until outside the galaxy. 

This assumption maximizes obscuration, since locking mass 
up in sub-resolution clumps would confine mass to smaller cover- 
ing fractions (see the discussion from simulations in Hopkins et al. 
|2005aT >. 



6.2 The Clumpy Torus 

In fact, we know that there must be sub-structure in the gas, be- 
cause cooling and star formation occur. Most of the mass in the 
ISM is probably locked into dense cold clumps. Unfortunately our 
simulation, limited by the physics included, does not predict the 
clump properties but only indirectly assumes an effective ISM state. 
However, with some simple assumptions, we can construct a sub- 
resolution estimate of all the relevant clump properties, without the 
introduction of any tunable parameters. 

Assume temporarily that most of the mass in the ISM is locked 
into N c \ dense clumps, with median mass M c i, size R c u and mean 
density p c \ = M c \/ (4tt/3)R 3 c1 . Define the density contrast p c i =xp, 
with respect to the volume- average background density p. We make 
two assumptions, both just at the order-of-magnitude level: that 
the clumps are quasi-virial, and that they are in pressure equilib- 
rium with the external medium. The first implies that whatever sup- 
ports the clump generates an effective pressure P c \ ~ Pc\Vc\ where 
Vc\ rsj GM C \/R C \. But this is justP c i = GT% h where E c i is the column 
density through the clump ~ p c \Rc\. To within a factor of two or so, 
this is even true for clumps in free-fall collapse, so is likely to be 
robust. We know the external effective (volume- average) pressure 
of the medium, P e ff - this is just the volume- average pressure used 
for all SPH calculations. It is straightforward to then set P c i ~ Peff, 
and obtain 



^cl — V Peff/ G . 



(19) 



Pressure equilibrium is a less certain assumption, but if we were 
to force a mass-radius or linewidth-radius relation similar to the 
observed Larson's laws in GMCs (a oc P 1 ^ 2 ), we would obtain 
the same dimensional scalings. 7 Assuming that clumps follow the 
Jeans mass and radius in a self-regulating Q—\ disk actually also 
results in the same dimensional scalings, so it may be robust in a 
variety of regimes. 

The probability of a path length Ar intersecting a cloud is 
given by p — (A/ci/Vtot) cr c \ Ar, where V to t is the total volume, and 
cr c i rsj ttR.Ii the clump cross section. But since N c \ ~ M tot /M c \, this 
simply reduces to 



p cl ~ pY>- x l Ar = p (Peff /G)~ 1/2 Ar 



(20) 



The only two quantities we ultimately care about, the proba- 
bility of intersecting a clump, and the clump column, have the use- 
ful feature that the clump density contrast and number of clumps 
completely cancel out. Thus, for any system where the mass is con- 
centrated in quasi-virial, pressure-equilibrium clumps, we can de- 
termine the column density distribution and probability of sight- 
lines seeing clumps based only on reference to well-determined 
volume- average gas properties in the simulations (p and P e ff). Of 



7 These clouds cannot, however, simply follow an extrapolation of the local 
GMC scalings. The local GMC size-mass relation implies an approximately 
constant clump surface density X ~ 10 22 cm -2 . But this is much less than 
the mean surface density of gas already at these radii, so any substructure 
must obey a relation at least different in normalization. 



course, the external pressure is set in part by our adjustable g eos , so 
it is important to examine the consequences of that choice. Higher- 
order detailed radiative transfer effects will depend on the specific 
clump sizes and other internal properties, but these are not our focus 
here. Because of the cancellation of the exact size and density con- 
trast (and correspondingly clump mass), the above relations hold 
for an arbitrary spectrum of clump masses, sizes, and/or densities. 

The column density along a given line-of- sight can then be in- 
tegrated outward from the BH. For each integration step Ar along 
the ray (taken to be increments of eh sm \, where e ~ 0.01 <C 1 and 
h sm \ is the local smoothing length at each point), we determine the 
probability p c \ of intersecting a clump, and probabilistically assign 
the ray a collision or not. If there is a collision, the integrated col- 
umn is increased by E c i. If not, the column is integrated through the 
"diffuse" (non-clump) phase of the ISM. The mass fraction in this 
phase (i.e. mass fraction not in star-forming clumps) is determined 
implicitly in the GADGET code (see Springel & Hernquist 2003), 
but is always small and should have near-unity volume filling fac- 
tor. 

Whether or not these assumptions are justified in detail, this 
provides a useful toy model, and we show that it can account for a 
number of observations. Moreover, on galactic scales, the assump- 
tions above have been borne out by a large number of independent 
observation s (|Larson|1981[|Ward-Thompson et al.|1994 [Sco ville 



|et al.|1987||Solomon et al.|1987||Ro solowsky p007l puller & My- 
|ers|19921|Andre et al.|1996||Blitz & Rosolowsky|2006| ). Of course, 
such clumps as observed locally could not survive the tidal forces 
near a supermassive BH. But even on nuclear starburst scales, it ap- 
pears that the star formation efficiency per clump dynamical time is 
low, implying they must be quasi-virial and not wildly out of pres- 
sure equilibrium ( |Tan et al.|2006l|Krumholz & Tan|2007t . Similar 
constraints come from clump structure in the narrow-line region 
(e.g. |Crenshaw en^|2000l |Rice et al.|[2006] >. And the fact that 
similar dimensional scalings arise from Jeans considerations im- 
plies they are likely to be generic to within factors of a few. Finally, 
we note that the dynamic range in column density is so large that 
violations of the above assumptions would have to be more than 
order-of-magnitude in order to qualitatively affect our conclusions. 

6.3 Column Densities: Model and Observations 

Figure [8] compares the resulting column density distribution, for 
simulations with varied q eos (each sampled at a random time near 
their peak of accretion). We compare the distribution from the 
"smooth" and "clumpy" torus models above, and that observed. 
Because of the dynamic range in Nu predicted, we are specifi- 
cally interested in comparison with samples sensitive to Compton- 
thick populations. We compile the (estimated intrinsic) distribution 
of column densities determined from the INTEGRAL/IBIS AGN 
sample of Ma lizia et al.|p009[ 20 — 40keV), the predominantly 
SWIFT/BAT sample o f|Treister et al.|p009l - lOOkeV), and the 
nearby OIII sample in Risaliti et al.| ( |199 9) (this is a Type 2-only 
sample, so we normalize to their estimated total fraction of Type 2 
AGN). The latter sample is most complete at the highest columns 
> 10 25 cm -2 ; none are sensitive to AGN with Nn > 10 26 cm -2 . At 
lower column densities, these are consistent with a wide variety of 
hard X-ray observations from e.g. Chandra and X MM (|Ueda et al.| 
|2003l[LaFranca et al.|2005[|Sirverman et al.|20051|Hasinger|2008] L 
And more recent, independent analysis of larger SWIFT/BAT sam- 
ples also agrees well ( [Burlon et al.|2011| ). 

Unsurprisingly, the predicted columns in the smooth torus 
model are uniformly large, in conflict with the observations. This 



© 0000 RAS, MNRAS 000, 000-000 



The Dynamical AGN Torus 1 3 




0.8 
0.6 
0.4 
0.2 



0.0 



Maliziaetal. 2009 - 
Triesteret al. 2009- 
Risaliti etal. 1999 " 



20 



21 



22 



23 



24 



25 



26 



27 



log( N H ) [cm- 



Figure 8. Top: Column density distribution predicted by our simulations. 
Each line is a simulation with different sub-grid equation of state (i.e. feed- 
back/pressure support in the gas), as in Figure [7] Here, we assume there is 
no sub-structure in the gas (i.e. gas is perfectly smooth below our resolution 
limits). The obscuration is clearly over-predicted. Middle: The same, but 
assuming the gas is clumpy. The clumps are assumed to be quasi-virial and 
in pressure equilibrium with the outside medium - this completely deter- 
mines the predicted Nn distribution with no free parameters. Bottom: Ob- 
servational estimates of the column density distribution, fr om|Malizia et al.| 
( p009l black), |Treister et al.|p009l red), and |Risaliti efaLHl999| blue). The 
cutoff in these samples at Nu > 10 26 cm -2 is a selection effect. Allowing 
for a simple clumpy gas model, without any tunable parameters, provides a 
good match to observed Nu distributions. Because of the effects of gravita- 
tional maintenance of h/R, and the global similarity of mass distributions 
shown in Figure^] there is relatively little dependence of the column density 
distribution of q eos • 



is not a problem of there being "too much" gas - recall that the 
actual total gas masses and gas densities predicted at these scales 
agreed well with those in observed AGN (Figure [TJ. What this 
shows is that it is not possible to reconcile the observed central 
masses, gas densities, and/or SFRs of AGN with their obscured 
fractions, without invoking some small-scale gas clumping. The 
problem cannot simply be that systems are observed at different 
states either - as pointed out in |Hicks et ah] p009| ), several ob- 



served optically un-obscured AGN have instantaneous near line- 
of- sight volume- averaged gas densities in < 1 — lOpc that should 
naively imply columns of Nh ~ 10 25_26 cm -2 , similar to our pre- 
dictions here without sub-resolution clumping. And indeed direct 
observations on this scale have argued for such clumping (Risaliti 
let al.|20021|Mason et al.|20061|Sanchez et al.|2006l|Nenkova et al.| 
|2008b[ |Ramos Almeida e t al. 2009] |Hoenig & Kishimoto||2009l 
|Deo etal.|2011) . 

The column density distribution predicted by the clumpy torus 
model, on the other hand, agrees well with that observed. The ba- 
sic features are easily understood: the small mass fraction in the 
diffuse ISM phase shifts the main peak in the Nh distribution to 
lower values. The tail towards larger Nh is caused by obscuration 
by clumps. The relative "flatness" of the tail is broadly expected for 
vertical profiles similar to those in Figure [3] 

Although the systems plotted differ in some subtle details, 
there is little dependence on the parameterization of stellar feed- 
back (our q eos parameter). Why should the column density distri- 
bution be so insensitive to stellar feedback? Most important are the 
factors discussed in § |4.3[ i.e. the contribution of gravitational heat- 
ing which keeps the disks somewhat puffed up, and means that the 
gaseous scale height does not scale as strongly with q eos as might 
otherwise be expected. 

There are also two handy 'conspiracies,' in the clumpy torus 
scenario, which make the predicted column density distribution pri- 
marily a function of global, rather than local parameters. In the 
(near-polar) regime where p c \ <C 1, it is quite difficult in any model 
to obtain a column density radically different from those shown. 
This is because, even if all the mass is locked in cold clumps, a 
column of at least ~ iq 21-22 cm -2 will arise just from diffuse, non 



star-forming galactic gas on much larger scales (see Hopkins et al. 
2005b 2006a). We do see some systematic difference in the lowest 
columns seen, because the exact mass in the "diffuse phase" de- 
pends on the sub-grid model - but for almost any reasonable model 
this mass is small, so these differences are all in the un-obscured 
range (and therefore dominated by or comparable to galaxy-scale 
effects). In the opposite (near-disk plane) regime, where p c \ > 1, 
the total column encountered is ~ E c i p c \ ~ p Ar - i.e. in the opti- 
cally thick regime the column density is simply the same as that of 
the "average medium," independent of the gas properties or phase 
structure so long as the global dynamical properties predicted are 
similar (physically, this simply represents where clouds will begin 
to overlap, thus making a more uniform molecular medium). This 
is true even if we discard our assumptions of virial and/or pressure 
equilibrium. It is only in the intermediate column regime (which 
interpolates broadly between the two, so we do not expect any fea- 
tures or particular sensitivity to appear) where the detailed assumed 
model of clump properties makes some difference. 

Figure [9] illustrates how the column density varies with incli- 
nation angle, (for the "clumpy" scenario). Qualitatively, the be- 
havior is expected: columns increase towards the disk plane. There 
is, however, significant scatter in the column density at a given 
even within a given simulation at a given time. Strikingly similar 
results are seen in simulations by Wada et al. ('2009), despite in- 
cluding a very different model for stellar feedback, and ignoring 
the role of self-gravity. We also show the expectation value of the 
number of clumps encountered along each sightline. As expected, 
this increases along the disk plane. Pole-on, it is ~ 0. 1 <C 1 in al- 
most all cases. Edge-on, it typically reaches ~ a few. 

These values are consistent with various indirect constraints 
from attempts to model AGN SEDs (M ason et al.|2006||Shi et al. 
|2006l|Tnompson et al.|2009l|Ramos Almeida et al.[2009[|Mor et al, 
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Figure 9. Obscuration properties from the clumpy-torus model in Fig- 
ure [8] versus viewing angle (cos# = is edge-on; cos# = ±1 is polar). 
For q eos survey. Top: Integrated column. Dense clusters of points at lower 
A^h pass through the diffuse ISM only, points scattered to higher Nu en- 
counter clumps along the line-of-sight. The appearance of bimodality is 
somewhat artificial (see Figure|8] for most simulations, there are not actu- 
ally two peaks). Each point is a sightline to the BH in a given simulation. 
Bottom: Expectation value for the number of clumps encountered along a 
given sightline. The clump number density is always Poisson, with ~a few 
clumps along the typical edge-on line of sight, and rapidly declines above a 
scale-height ofh/R~l / 2. 



2009 , Hoenig & Kishimoto 2009 ; Nenk ova et al.|2008a|b) . Almost 
universally, these studies have found that a similar clumpy torus is 
required, with a number of clumps of order several along the edge- 
on lines of sight, characteristic locations/outer radii of most of the 
clumps from ~ 1 — lOOpc from the BH, and (where constrained) 
radial clump distributions with roughly power-law scaling dp/dr~ 
We find, for our typical gas surface density profiles S gas oc 



-(0.5-1) 



i dp/dr 



-(0.7-1.2) 



over the dynamic range of interest 



here. 

The number of clumps can be crudely estimated from 
Eqn. [20| It is straightforward to show that this equation reduces 
to (Slumps) = / pddr ~ (R/h)Q~ 1/2 where h w c s /Q is the 



scale height of the torus and Q is the usual Toomre Q. For a 
self-regulating disk, therefore, with Q ~ 1, we naturally expect 
(Afdumps) ~ (h/R)~ l ~ a few. The same scaling pertains if we dis- 
card pressure equilibrium and instead assume clumps are character- 
istically Jeans-scale in a Q ~ 1 disk (since then the scale of clumps 
within R is ~ h). 

The characteristic value of a few clumps is also interesting be- 
cause it implies that one is almost always in the Poisson regime. 
This has several implications. First, there should be a large scat- 
ter between the column observed and actual viewing angle, con- 
sistent with a wide variety of observations (see references above). 
Second, clumping has a number of important radiative transfer ef- 
fects, which will be discussed in subsequent work. Third, this al- 
lows for highly variable obscuration. A clump moving through the 
line of sight can lead to variation in the column density by several 
orders of magnitude. The detailed variability will depend on the 
clump size spectrum and other properties, but the maximal vari- 
ability timescale should scale as ~ R c i/RQ(R); since most of these 
clumps are at ~ 0. 1 — 1 pc, the constraint that clumps not be tidally 
shredded (p c \ > Mbr/R 3 , and N c \ > 1) sets an upper limit to the 
variability timescale of ~ 5 — lOOyr (for 0.1 — lpc), for a 10 8 Mq 
BH. For partial obscuration, a more realistic clump density con- 
trast and/or larger clump number, the obscuration could vary on 
a timescale 0.01 —0.1 times this (i.e. months-year). Such rapid, 
extreme variability in X-ray obscuration has been seen in several 
AGN (|RisalitietaLf2 002 2005 ||Matt et al.|2003||Lamer et al.|2003| 



Guainazzi et al. 2005 ; Fruscione et al. 2005 ; Immler et al. 2003) 



7 THE OBSCURED FRACTION AND TORUS 
PROPERTIES 

We now examine how the column density distribution depends 
on global properties. For the sake of comparison with observa- 
tions, we parameterize the distribution by means of the "obscured 
fraction": specifically, the fraction above a given column density 
Nh > 10 22 cm -2 (a value typically adopted in observational stud- 
ies). Henceforth, we ignore the "smooth torus" model - it does not 
agree with observations and gives uninteresting (always near-unity) 
obscured fractions. 

Figure [l0| compares the obscured fraction in the "clumpy" 
model with a number of nuclear properties. For each simulation, 
we measure the relevant properties at randomly sampled times and 
viewing angles. We show the obscuration versus total mass inside 
some small radius (essentially, the BH mass), versus gas mass, ver- 
sus the nuclear SFR, and versus the BH Eddington ratio. 

Unsurprisingly, / bsc increases with the gas mass inside a 
small radius < lpc. Note, however, that the correlation is weak: 
/obsc oc Mg/ S 4 . The midplane columns should increase more rapidly 
with Mga S , but these are already optically thick - the obscured frac- 
tion grows slowly with the fraction of sightlines above the disk that 
(at higher column) become optically thick. More interesting is the 
correlation this implies - / bsc also increases with the nuclear SFR. 
Behavior along these lines has been observed at a wide variety of 
scales - Type 2 AGN are more likely to be found in more rapidly 
star-forming hosts, and/or hosts with younger stellar populations 
(|Brotherton et al.|1999||Canalizo & Stockton|2001||Yip et al.|2QQ4| 



Jahnke et al. 2004, Zakamska et al. 2006, Nandra et al. 2007 



2004, 



verman et al. 2008). The observational correlation appears to be 



particularly strong when the nuclear stellar populations are isolated 
<|Shi et al.|2 007 , Wa ng et aL|2(Xr7l|Imanishi|2002[|Imanishi & Wada| 
|2004| |Davies et al.||2007| ). Note that the SFRs inside of ~ lOpc 
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Figure 10. Obscured fraction (fraction of sightlines that encounter a column density Nu > 10 22 cm -2 ) for our simulations as a function of various properties. 
Each simulation is sampled at several random times both before, during, and after its period of peak activity, and we show results from all simulations 
with different initial conditions and q eos values. Points are colored by the instantaneous gas fraction inside of lOpc (from red at / gas < 0.05 to dark blue 
at /gas > 0.8). Top Left: Versus Mbh (we add ~ 0.1 dex scatter in Mbh so the points can be distinguished). For otherwise the same conditions (e.g. M gas ), 

— 1/2 1/4 

/obscured declines with Mbh as approximately oc M BU ' . Top Right: Versus M gas inside 1 pc, at fixed Mbh- Obscured fractions increase weakly (oc M gas ) with 
gas mass. Bottom Left: Versus the star formation rate inside lOpc. There is a weak increase in /obscured with the circum-nuclear star formation rate, driven by 
the dependence on M gas (slightly weaker oc 

^■15-0.2 

, because M* ' is super-linear in M gas ). Bottom Right: Versus inflow rate into small radii, at fixed 



A^bh = 3x 10 7 Mq (if this continued to arbitrarily small radii, this would be proportional to the Eddington ratio M^n/M^dd : 
is significant, but very weak (oc M 01 ). 



: L/L^dd)- The dependence here 



can reach large ~ IOMq yr values; however, as shown in Hopkins 
|& Quataert| ( |2010a| ) (Fig. 14), this is correlated with the BH in- 
flow rates on these scales as M* ~ M B h (again, both tracing the gas 
mass supply) - for a less extreme quasar the "zero point" expected 
inside ~ lOpc would be more like ~ 0. 1 — 1 Mq yr -1 . Conversely, 
ULIRGs and mergers with more pronounced star formation in their 
nuclei are more likely to host obscured Seyferts or quasars, whereas 
those with slightly older populations are more likely to exhibit Type 
1 sig natures (|Farrah et al.|2003] [20051 |Sanders|1999||Guyon et al.] 



[2006l|Dasyra et al.|2006[|Yuan et al.|2010] r 

Most likely, at least some of this trend owes to the role of 
AGN feedback in clearing away some of the gas and dust (see 



e.g.|Sanders et al.||1988a| |Hopkins et al.|2005c|b| [Hopk ins 201 1 



Gran ato et al.|2 004, Narayanan et al. 2006 ), but it can simply arise 
as we see here from the larger gas and dust supply "burying" the 
AGN until star formation exhausts much of that material. We stress 
that this is a true nuclear- scale (< lOpc) correlation here, and the 
nuclear SF contributes negligibly (< 0.3 %) to the total SFR; there 
is not necessarily any predicted correlation between the AGN ob- 
scuration and the total/large- scale galaxy SFR. 



Given similar gas properties, the obscured fraction decreases 
with BH mass. This is expected because the BH gravity provides 
a stabilizing force that tries to "flatten" the torus (for fixed gas 
properties, the disk h/R oc M^ 2 ). If more luminous AGN are, 
on average, more massive BHs, then this suggests an inverse cor- 
relation between QSO luminosity and obscured fractions. Indeed, 
the existence of such an apparent correlation is well-established 
([Hill et al.|1996l|Simpson et al.|1999[|Willott et al.|2000||Simpson 



& Rawlings 2000, Steffen et al. 2003, Ueda et al. 2003, Grimes 



et al. 2004 



Hasinger 2004, Sazonov & Revnivtsev 2004; Barger 
& Cowie|2005[|Simpson|2005[|Hao et al.|2005[|Gilli et al.|2007[ 
Hickox et al. 2007 , Hasinger 2008 ). However, it is still unclear pre- 



cisely how much of this correlation owes to alternative possibilities 
such as simple dilution by the host galaxy and/or differences in the 
Eddington ratio distribution and accretion state (for a detailed dis- 
cussion, see |Hopkins et a l. 2009b and references therein). More- 
over, without a full cosmological model to predict e.g. the distribu- 
tion of active BH masses and Eddington ratios, we cannot forward 
model the BH luminosity distribution to construct a direct compar- 
ison with observations. But the predicted scaling here is not es- 
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pecially strong; it may well be that additional physics is needed 
to recover the full observed correlation - most commonly, AGN 
feedback is invoked to "blow away" some of the torus in the most 
luminous systems (see references above). 



8 DISCUSSION 

We have studied AGN obscuration in a series of multi- scale hy- 
drodynamic simulations that can self-consistently follow gas from 
> lOkpc galactic scales to < 0.1 pc. These simulations include the 
full self-gravity of stars and gas (along with BHs and dark matter), 
gas cooling, and star formation, along with varied prescriptions for 
feedback from young stars; these are all critical to the behavior we 
see, and have not before been simultaneously modeled on nuclear 
scales. In these simulations, inflows from large scales, when suffi- 
ciently large, lead to a cascade of instabilities on small scales, ul- 
timately yielding large nuclear gas masses and accretion rates onto 
the AGN. The scenario is qualitatively similar to the "bars within 
bars" model, but there is a high degree of variability and morpho- 
logical diversity at each stage (with spirals, bars, clumps, flocculent 
structures, all present and alternatively powering inflows and out- 
flows) - a more apt description would be "stuff within stuff" ( |Hop-| 
|kins & Quataert 2010a). Once gas nears the radius of influence of 
the BH, it generically forms an unstable m—\ mode (a lopsided 
or eccentric disk) that slowly precesses about the BH (Hopkins & 
Quataert 2010b). The stellar and gas disk precess differently, lead- 
ing to strong gravitational torques that can drive accretion rates of 
up to rsj IOMq yr -1 onto the BH. 

In this paper, we show that these nuclear, lopsided disks in 
fact naturally account for the long-invoked "toroidal obscuring re- 
gion" used to explain the obscuration of Type 2 AGN. Up to now, 
these models have been essentially phenomenological - we show 
for the first time the formation of sub-pc scale obscuring structures 
from galaxy-scale inflows, and in a suite of ^ 100 simulations show 
that they are quite generic and arise ubiquitously with this inflow 
scenario. We show that the global dynamical properties - gas and 
stellar densities, density profiles, kinematics, gas fractions, and star 
formation rates, agree well with observations of AGN obscuring 
regions from scales as small as 0.1 pc to > lOOpc. 

This implies a fundamentally new paradigm in which to view 
the obscuring region or "torus." Far from being a passive bystander 
or simple fuel reservoir for the accretion process, it is itself the 
driver of that accretion. The torus is the gravitational structure 
on scales within the BH radius of influence that torques the gas 
and forces continuous gas inflow onto the BH. The same lopsided 
modes that drive accretion can also provide the scale height, col- 
umn density distribution, and characteristic gas properties of the 
structure. 

As such, the predicted torii have non-trivial substructure: both 
small-scale clumping in the gas (discussed below), global m—\ 
patterns, and warps/twists arising from bending modes at a range 
of radii. On large scales, the m = 1 modes tend to manifest as lop- 
sided/eccentric disks, or one-armed spirals; on small scales, they 
become more tightly wound spirals. Their typical amplitude in sur- 
face density is expected to be at the ~ 10' s of percent level (see 
Hopkins & Quataert 2010a); the amplitude of induced non-circular 
velocities and corresponding magnitude of "offsets" of the BH from 
the spatial center of the galaxy nuclei (in units of the BH radius of 
influence) are about the same. Maser observations may show in- 
dications of asymmetry in the structure around nuclei ( Schinnerer 
|et al. 120001 |Greenhill et al.|2003l|Kondratko et al.|2005l|Fruscione1 



et al. 2005, Kondratko et al. 2006a); it is also possible that mea- 
surements of the velocity structure of e.g. molecular emission lines 
from the torus region may be able to measure such asymmetries in 
the near future. 

We have argued that a large number of obscuration properties 
traditionally associated with "feedback" processes from AGN and 
star formation may, in fact, be explained by purely gravitational 
physics. Even in the absence of feedback, properly including the 
full self-gravity of gas and stars leads to disks with large h/R, suffi- 
cient to account for the observed column density distribution. This 
arises because of a combination of large-scale warps and twists (for 
example, where the lopsided disk mode meets an outer bar) and 
bending modes within the disk itself. The latter will, even in the 
absence of any large-scale twists or warps, tend to pump up h/R 
wherever the eccentric disk mode is excited until an order unity 
h/R ~ \a\ is reached. Since bending modes are fast modes (pattern 
speed ~ Q), this can continuously transfer energy from the orbital 
motion to vertical motions on a single dynamical time, maintain- 
ing vertical scale heights even when the cooling time is arbitrarily 
short. 

These warps and twists also naturally lead to the observed lack 
of correlation between nuclear- scale disk inclination angles and 
those of their parent/host galaxies. This will be even more promi- 
nent in systems which are driven on large scales by mergers, but can 
occur even in entirely secularly fueled AGN. They also account for 
observed gas velocity dispersions in AGN nuclei, and the correla- 
tions between those dispersions and quantities such as the local gas 
mass, star formation rate, and mass in young stars (all via their in- 
herent dynamical correlations, not via any feedback channel). The 
efficiency of gravitational torques and induced inflow also natu- 
rally leads to convergence in nuclear gas masses and density pro- 
files, leaving relic "cusps" similar to those observed (Hopkins & 
|Quataert|20lTb) . 

These mechanisms can naturally explain observed global 
quantities such as the gas scale heights, masses, and density pro- 
files. However, modeling the actual obscured fraction of AGN re- 
quires a more explicit model for the actual sub- structure on Jeans 
mass scales and well below, in the ISM surrounding black holes. 
We show in fact that any model which matches the observed dy- 
namical properties (particularly global gas masses), but assumes 
"smooth" gas (uniformly distributed, say, out to some scale height 
corresponding to the average obscured fraction), will simultane- 
ously fail to explain the observed column density distributions. A 
natural explanation for this discrepancy is that the gas is clumpy on 
multiple scales, broadening the column density distribution along 
all sightlines. There must, in fact, be structure on the relevant 
scales, since we know there is star formation at these radii (so some 
gas must reside in dense, tidally bound star-forming clumps). Un- 
fortunately, our present models do not explicitly resolve the nec- 
essary physics of star formation and GMC formation/destruction 
via stellar feedback needed in order to explicitly simulate the sub- 
structure of the gas down to these scales. 

However, we find that we can obtain predicted column density 
distributions in good agreement with those observed if we assume 
that whatever sub-grid clumps exist obey a couple of basic assump- 
tions: namely that they are (at least at the order-of-magnitude level) 
near both virial and pressure equilibrium (or, instead of pressure 
equilibrium, that they are Jeans- scale in a Q ~ 1 self-regulating 
disk). These assumptions are sufficient to (statistically) predict the 
column density distribution that would be observed, regardless of 
the actual clump mass spectrum and physical origin (and without 
any adjustable parameter introduced). The predictions agree quite 
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well with the column density distribution of both un-obscured, ob- 
scured, and Compton-thick AGN. This suggests that these basic 
properties should still hold for substructure in these regions, and 
that - if so - the uncertain quantities in our simulations (such as 
the feedback prescription and star formation recipe), make no dra- 
matic difference in the column density distribution, since its key 
properties are set by the basic dynamics above. Essentially, if these 
assumptions hold, the distribution of observed column densities to- 
wards AGN is itself a natural consequence of gravitational clump- 
ing at the Jeans length/mass in a self-gravitating, globally Q ~ 1 
disk - no exotic wind physics (driven by either stars or AGN) need 
to necessarily be invoked. 

Higher-order probes of the structure in this region, for exam- 
ple studying the clump properties (their sizes and masses), con- 
straining the ratio of stellar feedback to dynamical support in driv- 
ing scale heights, and making predictions for line structure and 
other effects that might be used to probe the sub- structure and lop- 
sided precession that power accretion, will require detailed treat- 
ment of the radiative transfer from the accretion disk through the 
circumnuclear region. This will be the subject of a future paper, 
and should enable a host of new predictions for comparison with 
future observations. 

Another important next step will be the inclusion of realistic, 
physically motivated feedback models. Coupling our simulations 
with radiative transfer will be a major advance. Although this ap- 
proach will not be strictly self-consistent, we will, for the first time, 
be able to examine how radiation pressure impacts inflowing and 
star-forming gas using a realistic description of multi-scale AGN 
gas distributions from ^0.1 — lOOOpc scales. In particular, to study 
where the photon momentum is absorbed (compare Murray et al. 
|2005[|Ciotti et al.|2010l|Hopkins et al.|2011a| ), how radiation pres- 
sure profiles vary throughout the gas, how photon diffusion may 
affect the role of feedback ( Thompson et al. 2005 ), and whether a 
realistic clumpy gas medium suppresses or enhances the efficiency 
of feedback-induced "shutdown" in star formation (Hall et al. 2007 ; 
Hopk ins & Elvis|2010|[Tortora et al.|2009| >. In bright quasars and/or 
nuclear starbursts, the gas structure may well be modified not just 
locally (in terms of its dumpiness or sub-grid pressure support), 
but globally by strong outflows driven, for example, by radiation 
pressure (see Debu hr et al.|201l[ and references therein). Even in 
the regime where some material is being expelled at the escape 
velocity, it is difficult to alter many of the basic dynamical proper- 
ties of the gas (total mass enclosed and its relation to inflow rates, 
obscured fractions, etc) at the order-of-magnitude level (see e.g. 
Marc oni et al.|2008] >, but may well make a large contribution to the 
observed scale height of obscuring material and can be critical to 
understanding how AGN self-regulate, why torii exhibit complex 
sub- structure, and perhaps scalings of obscured fraction with lumi- 
nosity and/or redshift. 

We have focused here on small-scale obscuration, at radii tra- 
ditionally associated with the AGN "torus." We stress, however, 
that this does not mean that there is a single object that accounts for 
the obscuration of all systems. As is evident in all of our compar- 
isons, the gas distribution is truly continuous. Of course, there will 
be gas on small scales near the BH whenever it is active, which can 
occult and obscure different emission regions. This may take the 
form of an AGN wind, especially in high Eddington ratio systems 
(see e.g. |Elvis|2000l|Elitzur & Shlosman|2006] ). 

There is also well-resolved gas from the host galaxies in these 
systems. The latter, on say > lOOpc scales, is not likely to be 
Compton-thick, simply because the characteristic Jeans scales, etc. 
are too large. However, this can easily dominate the production 



of more moderate column densities Mi ~ 10 cm - . This "host 
galaxy obscuration" is especially important in the early phases of 
inflow forming a kpc-scale starburst, as the central kpc may be 
isotropically enshrouded in dust for a time ~ 10 8 yr. Such obscu- 
ration as it arises in simulations of galaxy mergers has been dis- 
cussed at length in e.g. Hopkins et al. (2005b 2006a); Hopkins 
|& Hernquis"t| ( |2006| ); [Hay ward et al.| ( |2011| >, and we refer to those 
papers for more details. Observations have also made it clear that 
a significant fraction of obscuration must come from host galax- 



ies (especially in starbursts and edge-on disks; see e.g. Zakamska 



et al.] 2006 ; Rig by et al.[20 06 ; Martinez-Sansi gre et al.|20091|Lagos 
et al. 2011) Since our current simulations are the first to simulta- 



neously resolve both the nuclear scales where very large columns 
arise, and the galaxy scales where more moderate but potentially 
more isotropic (or at least differently-oriented) columns can arise, 
it will be interesting to investigate the relative contributions to ob- 
scuration from different scales, as a function of evolutionary stage 
and galaxy/BH properties. Because the AGN spectrum changes as 
it moves through the inner obscuring regions, this will require a full 
treatment of radiative transfer, as described above. 
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